Nonlinear Scattering of a Bose-Einstein Condensate on a Rectangular Barrier 
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We consider the nonlinear scattering and transmission of an atom laser, or Bose-Einstein con- 
densate (BEC) on a finite rectangular potential barrier. The nonlinearity inherent in this problem 
leads to several new physical features beyond the well-known picture from single-particle quantum 
mechanics. We find numerical evidence for a denumerably infinite string of bifurcations in the trans- 
mission resonances as a function of nonlinearity and chemical potential, when the potential barrier 
is wide compared to the wavelength of oscillations in the condensate. Near the bifurcations, we 
observe extended regions of near-perfect resonance, in which the barrier is effectively invisible to 
the BEC. Unlike in the linear case, it is mainly the barrier width, not the height, that controls the 
transmission behavior. We show that the potential barrier can be used to create and localize a dark 
soliton or dark soliton train from a phonon-like standing wave. 

PACS numbers: 03.75.-b, 03.75.Pp, 03.75.Lm 



I. INTRODUCTION 

Bose-Einstein condensates (BECs) in the dilute gas 
limit and in the presence of many kinds of external poten- 
tials are modeled very well by the nonlinear Schrodinger 
equation (NLSE) [1, 2]. BECs with potential barri- 
ers have many practical applications, including atom 
lasers [3-5] and high-resolution holography [6, 7]. Two 
recent experiments at Washington State [8] and Rice [9] 
have considered the dynamics of a finite barrier dragged 
over an effectively ID BEC, giving one the opportunity to 
investigate the Landau criterion [1] under the constraint 
of a ID system. Instead of the quantized vortices com- 
mon to supcrfluids in two and higher dimensions, dark 
solitons appear in ID. As the stationary NLSE with a 
delta function potential can be solved exactly [10-13], it 
is convenient for theory to model localized potential bar- 
riers with delta functions. For instance, multiple delta 
functions have been used to investigate intriguing prob- 
lems such as nonlinear band structure [14] and Anderson 
localization [15]. However, as we show in this article, 
when the delta function is replaced with a potential of 
finite width new features arise, including a denumerably 
infinite series of bifurcations in the transmission reso- 
nances, and one or more dark solitons localized on the 
barrier. 

Since the stationary NLSE can be solved exactly for 
any piece-wise constant potential, we take the potential 
of finite width to be rectangular in form [16]. In single- 
particle quantum mechanics, stationary solutions of the 
Schrodinger equation for a rectangular potential barrier 
or well is a classic problem taught in undergraduate quan- 
tum mechanics courses. The nonlinearity inherent in the 
mean-field description of the BEC as captured by the 
NLSE substantially changes this well-known problem, as 
we will show. The physical context of our problem is 
the steady-state transmission behavior of an atom laser 



incident on a barrier. We assume that the BEC is con- 
fined in the transverse directions by a harmonic oscillator 
trap, and that its behavior is quasi-one-dimensional [17- 
20]; that is, the longitudinal direction of the BEC is 
much larger than the transverse directions and the heal- 
ing length, and the chemical potential is much larger 
than the transverse excitation energy. Further, we as- 
sume that the width of the potential barrier is much less 
than the longitudinal dimension of the BEC, so that the 
system is effectively longitudinally infinite and far-field 
effects may be neglected. This physical situation was ex- 
perimentally produced in the Washington State and Rice 
experiments [8, 9]. 

Past studies of the ID stationary NLSE with vari- 
ous boundary conditions and potentials are very exten- 
sive, including both repulsive and attractive interactions; 
see [21] and references therein. In this article we focus on 
the more common case of repulsive interactions, although 
our techniques work for general interactions in the NLSE. 
Our own work has often treated cases described by piece- 
wise constant potentials, starting with a uniform poten- 
tial with periodic or box boundary conditions [19, 22, 23]. 
This work was later extended to a discontinuous poten- 
tial step and to a delta-function potential [12]. Other 
studies in this direction included the Kronig-Penney po- 
tential [14] and the bichromatic lattice [24]. 

The potential barrier or well has been treated previ- 
ously in restricted circumstances in three previous stud- 
ies. First, Carr, Mahmud, and Reinhardt [25] found par- 
ticular classes of bound states in the potential well. Such 
bound states are localized or partially localized, and it 
was found that experimental parameters could be tuned 
to achieve different regimes of tunneling. Second, Ra- 
pedius, Witthaut and Korsch [26] took the region outside 
the well to be strictly linear. They found a bistable be- 
havior in the transmitted flux for unbound states. The 
critical nonlinearity at which bound states appear was 
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discovered. Third, Ishkhanyan and Krainov [27] consid- 
ered the limit of very small nonlinearity. They used per- 
turbation theory to expand the NLSE, and a multi-scale 
method to find the solutions. The reflection coefficient R 
could then be determined by the usual methods of linear 
quantum mechanics, decomposing the wave functions in 
each region into left- and right-traveling waves via the 
superposition principle. The latter two studies had in 
common that they retained the concept of the superposi- 
tion principle outside the well. We discard this concept in 
favor of full nonlinearity, solving the complete nonlinear 
problem without approximations of any kind. Additional 
related work on scattering of a bright soliton on a poten- 
tial barrier or well has treated a variety of applied and 
foundational issues in quantum mechanics, from bound 
state spectroscopy and resonant trapping [28] to macro- 
scopic superpositions [29, 30] (Schrodinger cat states). 



In all studies of the stationary ID NLSE solitons play 
a key role. In the Washington State experiment [8], the 
barrier was produced by an elliptical laser beam, then 
dragged through the BEC. For intermediate drag speeds 
a train of dark solitons was observed in the presence of 
the barrier. Such solitons have been produced by other 
experimental methods, for example, by merging two co- 
herent BECs [31]. An experiment by Weller et al. at Hei- 
delberg studied the dynamical behavior and interactions 
of such solitons [32]. Solitons are localized, persistent, 
robust nonlinear structures which appear often in BEC 
experiments [33, 34]. A key feature of solitons in quasi- 
1D in the absence of an external potential is that they 
interact elastically and do not dissipate [35, 36]. Such be- 
havior was observed in the Heidelberg experiment, and 
was found to agree with numerical simulations of NLSE 
solution dynamics [32]. 



Our article is outlined as follows. In Sec. II we 
present the fundamental equations and a general solu- 
tion method. We emphasize that our work is completely 
analytical and symbolic up to one numerical integration 
used to characterize transmission; the NLSE solution it- 
self is analytical and is determined rigorously [37]. The 
appendices contains some key brief proofs to support this 
rigor, as well as numerical studies to support our analy- 
sis and a brief overview of Jacobi elliptic functional rela- 
tionships used in our analysis. In Sec. Ill we apply our 
solution method to scattering on the rectangular barrier 
and obtain characteristic solutions, highlighting cases for 
which solitons are localized on the barrier. In Sec. IV we 
study transmission of the BEC for both wide and nar- 
row barriers. In Sec. V we focus in on the transmission 
resonances, where the barrier is invisible, and discover 
a series of bifurcations in such resonances. Finally, in 
Sec. VI we conclude. 



II. EQUATION AND METHODS 
A. The ID Nonlinear Schrodinger Equation 



The ID NLSE is 



= iftJU(i,t), (i) 



where g = 2a s huj is the interaction strength or nonlin- 
earity renormalizcd to ID [19], a s is the s-wave scatter- 
ing length for binary contact interactions, co is the os- 
cillation frequency of the transverse trap as described in 
Sec. I, M is the mass of the atoms or molecules that are 
Bose-condensed, and V{x) is an external potential. We 
take the Landau interpretation of the wavefunction or 
order parameter: p = |\1/| 2 is the local BEC number den- 
sity and v = (h/m)dxAvg(^/) is its local velocity. Here 
and throughout this paper, a tilde denotes a dimensional 
quantity. 

Non-dimensionalizing (1) by scaling everything to har- 
monic oscillator units leads to the dimensionless or scaled 
NLSE, 



\^+9\Hx,t)\ 2 + V{x) 



where 



x = 



ojt, 



9 = 



# = *£ 1/2 , 



(3) 
(4) 
(5) 

(6) 
(7) 



and t = yjh/(muj) is the harmonic oscillator length. 

Equation (9) can be solved [22] for stationary states of 
the form 



^(x,t) = \fp{x)e 



\[4>(x)-nt] 



(8) 



where /i is the eigenvalue of the stationary ID NLSE. 
Then the ID NLSE becomes 



-\^+9\Hx,t)\ 2 + V{x) 



V(x,t) = ^(x,t). (9) 



We work with the non-dimensionalized stationary ID 
NLSE, Eq. (9), throughout the rest of our article. 
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Exact Solution for a Constant Potential 



Calculation of Transmission 



When the potential V(x) = Vo is a constant, one can 
prove that the dimensionless density p = 2 and the 
phase (j)(x) have the form 



We consider a potential barrier of form 



p(x) = Asn (bx + 6 \m) + B, 
dcj) a 
dx p ' 



(10) 
(11) 



where A is the density scaling, b is the translational scal- 
ing, m is the elliptic parameter, So and B are offsets, and 
a is an integration constant. The function sn is one of 
twelve Jacobi elliptic functions [38], which can be inter- 
preted geometrically as the elliptical analog of the cir- 
cular trigonometric functions. In this interpretation, the 
square root of the elliptic parameter m represents the 
eccentricity of an ellipse. In the limit that m — > 0, the 
Jacobi functions become the circular functions; in the 
limit m — » 1, the Jacobi functions become the hyperbolic 
functions. 

Relationships between parameters may be obtained by 
substituting the solutions (10) and (11) into Eq. (9) and 
equating coefficients of linearly independent powers of 
Jacobi elliptic functions. The relevant proofs of linear 
independence are given in Appendix B. We obtain the 
relationships 



A 

b^ 

^[b 2 + (A + 3B)g] + V a , 



a z = B(A + B)(b 2 + Bg). 



(12) 

(13) 
(14) 



Since a appears only in its square, and p(x) does not 
depend on a, all solutions with a =/= arc doubly de- 
generate. That is, ±a result in solutions with the same 
density p(x) and eigenvalue p, but phases of opposite 
sign. 

Mathematical and physical considerations lead us to 
conclude that not all possible parameter values are rele- 
vant to this problem [37]. The relevant parameter space 
is summarized in Table I. 



TABLE I: Mathematical conditions on both NLSE pa- 
rameters and solution parameters. 



Parameter 


Condition 


A 


Real 


B 


Real 


b 


Either Re(6) = or Im(6) = 


So 


No constraints 


m 


Real 


9 


Real 


M 


Real 


a 


Real 



V(x) 



0, x < x\, 

Vo, x\ < x < x 2 , 

0, X 2 < X. 



(15) 



We will take V) to be positive definite. The barrier can 
have any width and height, with the boundaries x\ and X2 
located at any positions along the x-axis. We define re- 
gions I and III to be the left and right sides of the barrier, 
respectively, and region II as the region over the barrier. 
We apply our fully general solution from Sec. II B to a 
numerical study of the transmission of the BEC across 
the barrier. We note that the solution and our code can 
be generalized to arbitrary piecewise-constant potentials 
with a finite number of jump discontinuities. 

Using the solution to the constant potential one can 
treat any piecewise constant potential by the use of 
appropriate boundary conditions, similar to the well- 
known method for finding stationary solutions of the lin- 
ear Schrodinger equation with a rectangular well/barrier. 
Our explicit boundary conditions, expressed in terms of 
density and phase, are 



dp(x 



dx 



dp(x) 



dx 



(j>(xf) + 27m, neZ, 



a 1 



(16) 
(17) 

(18) 
(19) 
(20) 



where Xi is the location of the zth boundary and ± indi- 
cates the right /left side of the boundary [39] . 

The effective potential experienced by the BEC in the 
region over the barrier is 



V eS = V +gp{x). 



(21) 



In the regime p > max(14fj), transmission over the bar- 
rier is classically allowed. For p < max(T4ff), transmis- 
sion is classically forbidden. In both regimes, our scat- 
tering displays quantum or wave-like effects, modified by 
the nonlinearity. The effective potential is why the lin- 
ear and nonlinear problems are so different, in addition 
to the lack of a superposition principle. 

We briefly remind the reader of the linear case. For g = 
only in Eq. (9), the wave function in each region can be 
split into a sum of two terms representing left- and right- 
traveling waves, using the principle of superposition. We 
can take the left-hand side as the incident one. On this 
side, the solution contains both right- and left-traveling 
waves, i.e., incident and reflected waves, with only right- 
traveling waves on the transmission side of the barrier. 
In this case, we may define the transmission coefficient 
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D. Numerical Methods 



T = 



(l^transl 2 ) 
(l*inc| 2 ) ' 



(22) 



where \I/trans is the transmitted wave function and 
is the incident wave function. The angle brackets, ( • ), 
denote an average value over one period of the func- 
tion. The definition given by Eq. (22) is standard in 
linear quantum mechanics. In this interpretation, T < 1 
over the entire domain of the system, and T represents 
the probability that a given particle will be transmitted 
across the barrier. 

However, in the nonlinear case superposition does not 
apply. We cannot define separate left- and right-traveling 
waves in this case, and thus the transmission coefficient 
is defined simply as 



(I* 



(l*/| 2 > 



{Pin) 
(Pi) 



(23) 



where and ^ m are the total wave functions in re- 
gions I and III, with pi and pm the corresponding num- 
ber densities. In Eq. (22) this would be equivalent to 
replacing the denominator with the sum of incident and 
reflected amplitudes and then squaring to get the total 
probability density to the left hand side of the barrier. In 
the nonlinear case, since neither region I nor region III 
can be said to be the incident side, we could just as easily 
have replaced the definition in Eq. (23) with its inverse. 

Thus, in the nonlinear case, we find that T may exceed 
unity, according to Eq. (23). Physically, output cannot 
exceed input, and we conclude that the "transmission" 
coefficient as we define it contains information for atom 
lasers incident on cither side of the barrier. Since we 
cannot choose only right- or left-traveling waves in the 
solution, due to the nonlinearity, we cannot define T to 
restrict incidence to only one side of the barrier, and must 
consider both circumstances in the same solution set. 

To calculate the average over the densities, we note 
that the period of sn 2 (bx + So\m) is 2K(m)/b, where 
K{m) is the complete elliptic integral of the first 
kind [38]. Thus we obtain the average density by 



(p) 



2K{m) 



J dx p{ 



x), 



(24) 



where the integral is taken over one period of p{x). Using 
the properties of Jacobi functions and elliptic integrals, 
one can show that if < m < I, 



(P) =B + A 



1 



E(m) 



m mK(m) 



(25) 



where E(m) is the complete elliptic integral of the second 
kind. If m ^ [0, 1], we must first apply a transformation 
to write the density in terms of Jacobi elliptic functions 
which depend on a parameter m' € [0, 1], then use the 
methods above to simplify the integral. The relevant 
transformations are given in Appendix A. 



We use Mathematica to compute the transmission co- 
efficient, Eq. (23). We require an internal precision of 100 
digits; our main reason for using Mathematica is its fea- 
ture of arbitrary internal precision, which is absolutely 
necessary for the highly singular problem of nonlinear 
scattering. Besides internal precision, we use a uniform 
numerical tolerance of I0 -5 for parameters. If a quan- 
tity is smaller than this value, it is taken to be zero; for 
example, we take g < 10 to be zero. We consider only 
repulsive interactions, i.e., g > 0, in our numerical anal- 
ysis. In addition, we restrict analysis to positive barriers 
Vb > 0, for consistency with the idea of scattering by an 
atom laser. However, as mentioned previously, our math- 
ematics are completely general, and our solution may be 
applied to attractive interactions and potential wells. 

To find a scattering solution, we take A, B, and So 
on the left side of the barrier, denoted by a subscript /. 
These parameters can be calculated uniquely from the 
average number density, momentum, and energy for a 
physical atom laser, as described in [12]. We consider 
the physical quantities g and p, as well as the barrier 
parameters, to be input parameters for an experiment, 
and take them as known. Other parameters on the left 
side of the barrier may be obtained from Eqs. (1 2)-(f 4). 
We then use the physical boundary conditions (20) to 
solve for parameters in regions II and III in terms of the 
known input parameters from region I. 

Our code is completely symbolic, with the exception of 
one numerical integration used in computing the trans- 
mission coefficient. In order to keep the code symbolic 
wherever possible, we used Mathematica's pure function 
routines. These are functions whose arguments are de- 
fined in terms of their position rather than being given a 
specific variable name. This construction avoids the pos- 
sibility of errors arising from multiply-defined variables, 
while still allowing us to maintain a consistent convention 
for functional dependencies. 

The analytical solution process using boundary condi- 
tions yields twelve solutions for (A, B) in regions II and 
III. Six of these are extraneous solutions, obtained as a 
result of squaring both sides of an equation. Since they 
are not true solutions, they do not satisfy the bound- 
ary conditions, and this is used as a filter in the code to 
discard these extraneous solutions. The remaining six so- 
lutions all correspond to the same functional form of the 
density, and we need only consider one. In our construc- 
tion of Eq. (10) as the polar form of a complex function, 
we may take p(x) to be real without loss of generality. 
Solutions with Im(p) ^ are mathematically invalid, and 
these are encountered in the code as a result of one or 
more assumptions breaking down. Therefore, we take 
only real solutions for (A, B). 

To obtain the average densities in Eq. (23), we numer- 
ically integrate the densities using Simpson's rule [40]. 
The transmission coefficient is computed as in (23), and 
is treated as a resonance if it is within an interval of 



±10 -4 around unity. 

In the linear limit g — > 0, exact analytical expressions 
may be obtained for the boundary conditions. We use 
these exact expressions, rather than Mathematica's Solve 
routine, in the limit of small g. By doing this, we avoid 
problems with infinities due to internal Mathematica pro- 
cesses. 

We also made use of several other methods to ensure 
the correctness of the computations performed by Math- 
ematica. In our experience, Mathematica does not al- 
ways handle \J — 1 correctly. Therefore, we kept j as a 
symbol, using replacement tables to handle powers of 
i, and avoiding Mathematica's internal complex- number 
routines wherever possible. In addition, we included sev- 
eral identities for Jacobi elliptic functions via replace- 
ment tables, as Mathematica's default processes do not 
make use of these identities for simplification. Finally, 
rather than relying on Mathematica to handle the special 
cases of complex argument and m > 1 internally, we used 
known transformations to write equivalent expressions in 
terms of Jacobi functions with real argument and param- 
eter smaller than unity. These expressions were used in 
the code for computations. 

Convergence was verified using several methods. These 
are detailed in Appendix C. 



III. SOLITON LOCALIZATION 

A large variety of solution types appear after following 
the methods of Sec. II. We present here some particu- 
larly interesting cases which are very far from the kinds 
of solutions found for the linear Schrodinger equation. 
We recall that, in the linear Schrodinger equation, the 
wavelength of the plane wave components across regions 
I, II, and III is the same; it is only the spatial offset in the 
arguments of the exponentials and the amplitudes that 
can be different. In the nonlinear case the wavelength 
can be quite different in all three regions. We present 
two such cases, in which a phonon-like standing wave is 
transformed into a strongly localized soliton or soliton 
train over the barrier, and then returns to a phonon-like 
profile. We recall that phonons are a limit of Bogoliubov 
quasi-particles for BECs, and appear as a small modu- 
lation on top of a large offset, qualitatively speaking, as 
ripples on a pond. Since we work with the fully nonlin- 
ear problem our phonons are not constrained to be small 
oscillations, but continue to appear as ripples of ampli- 
tude A and translational offset So on top of an offset of 
magnitude B. 

Figure 1 shows the density and phase for a local- 
ized soliton. A well-localized solution can be seen above 
the potential barrier. The barrier, indicated with a red 
dashed box in the figure, has height Vq — 1 and width 
20. The nonlinearity is g — 2.02 and the chemical po- 
tential fi = 2.404. We choose Ai = 1, Bi = I, and 
(5 / = 0. In regions I and III, outside the barrier, the 
elliptic parameter m< 1. In region II, over the barrier, 




FIG. 1: (Color online) Density (upper panel) and phase 
(lower panel) of a solution to the NLSE in which an 
incident near-sinusoidal wave produces a single dark 
soliton localized over the barrier. The three colors in the 
phase curve denote the three scattering regions. The red 
dashed curve shows the potential barrier for reference. 



m approaches unity and the peaks of the sn 2 function 
broaden significantly. The local minimum of the density 
over the barrier is a dark soliton, called "grey" because 
it does not have a node. Such a soliton is always mov- 
ing. Examining the phase profile, we observe two areas 
in region II: away from the dark soliton is a background 
superflow with a shallower slope, while over the soliton 
there is a characteristic phase jump displaying a sharp 
increase in slope. In this solution the superflow flows to 
the right with velocity v = d x <p(x), while the dark soliton 
moves to the left with an equal and opposite velocity; the 
two velocities cancel each other, leading to a stationary 
state. 

The phase has an overall linear envelope on either side 
of the barrier, with clearly visible oscillations on this 
background. The background slopes differ slightly be- 
tween regions I and III, and the slope over the barrier 
in region II is smaller by a factor of 3. This shows that 
the velocity profile of the BEC need not be the same on 
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pipe) 




20 40 60 

(a) 

<P(x) 




20 40 60 80 

(b) 

FIG. 2: (Color online) Density (upper panel) and phase 
(lower panel), similar to Fig. 1, but with a soliton train 
localized on the barrier instead of a single soliton. The 

input parameters are identical to those of Fig. 1; only 
the barrier width has been changed, from 20 to 90. We 

observe a train of 9 very deep solitons localized on the 
barrier. 



either side of the barrier. Although we fixed parameters 
on the left-hand side in order to find this solution, as 
described in Sec. II, from the physical perspective of an 
atom laser the BEC is incident on the right-hand side. 
Thus, as happens approximately half the time in our solu- 
tion method, we have in fact fixed the output rather than 
the input parameters. The overall transmission right to 
left is about 50%, and T > 1 according to our definition 
from Eq. (23). 

With a longer barrier of width x% — X\ = 90 and the 
same nonlinearity g = 2.02, chemical potential p — 2.404, 
and input parameters Ai = 1, Bj = 1, and Sqi = fixed 
on the left hand side, we obtain a soliton train localized 
on the barrier, as shown in Fig. 2. In Fig. 2 we again 
display the potential with a dashed red line for reference. 
The location of the dark solitons can be observed in both 
the density and the phase, similar to Fig. 1. Multiple 
solitons have been observed in a variety of BEC exper- 
iments [8, 32, 41]. However, such solitons have so far 



not been phase-locked into a train in BECs. Our barrier 
method constitutes a totally new way of producing multi- 
ple solitons locked into a train. The length of the barrier 
can be used to control the number of solitons in the train. 
We note that dark soliton trains have been produced in 
nonlinear spin waves in thin magnetic films [42, 43], but 
not via our scattering technique; also spin waves are in 
fact a damped driven system which are best modeled by 
an open-system version of the NLSE quite different from 
Eq. (9) [44]. 



IV. ATOM LASER TRANSMISSION 

We turn now to a more systematic exploration of the 
solution space. The nonlinear Schrodinger equation has 
a substantially larger parameter space to explore than 
the linear case. We divide this exploration into wide and 
narrow barriers. The barrier size xi — x\ can be compared 
to the BEC healing length, £ = 1/ y / 8Tr(p)a s , with (p) 
the average linear number density. However, the healing 
length only provides a useful comparison for low-energy 
excitations, as it describes how a uniform ground-state 
BEC is perturbed by a localized potential, e.g. a hard 
wall. This corresponds to well-separated dark solitons, as 
in region II in Figs. 1 and 2. As we treat a wide variety 
of excitations, many of which take standing-wave form 
which is very far from the ground state, e.g. regions I and 
III in Figs. 1 and 2, the correct comparison to identify 
wide and narrow barriers is in fact the wavelength of the 
excitations, 

= 2K(m) = 2K (2( M -y ) A (A+3B) 8 ) 
b ~ ^2(p-Vo)~(A + 3B)g' 

as can be calculated from Eqs. (12) and (10). Thus X2 — 
x\ <C Xj is the narrow barrier case and X2 ~ x\ ^> Xj is 
the wide barrier case, where the subscript j £ {I, II, III} 
refers to the region. For sufficiently small A all barriers 
are effectively wide. Holding other parameters fixed, in- 
cluding the nonlinearity, the wide barrier limit occurs 
for large chemical potential, since the wavelength is a de- 
creasing function of /i, as can be verified by consideration 
ofEq. (26). 

We found that the barrier height Vq is less important, a 
statement which we will support further in Sec. V. This 
dependence on width more than height is another sense 
in which the nonlinear case is very different from the lin- 
ear case; in the latter it is only the effective area of the 
barrier, 2MV (x2 — Xi) 2 /h 2 , that is important. The de- 
pendence on width arises mainly from the fact that the 
density in the nonlinear case can have different wave- 
lengths in different regions for the same solution, in con- 
trast to the linear case, as previously mentioned. 

For simplicity we will keep the input parameters Aj, 
Bj, and 5 j fixed to the same values as those used in 
Figs. 1-2 and Sec. Ill; after a massive exploration of the 
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FIG. 3: (Color online) Wide barrier, small nonlinearity. 
Transmission for wide barrier with increasing 
nonlinearity from bottom to top, starting with the 
linear case, g — 0, in steps of Ag = 0.01. Even a very 
small nonlinearity changes the problem because the 
left-right symmetry of the effective potential is broken. 



parameter space, we found all results to be qualitatively 
similar to those described below. 



A. Wide Barrier Case 

We consider a barrier whose width is much larger than 
its height. We take a barrier of width x-i — x\ = 20 and 
height Vo = 1 for the purposes of illustration. In Figs. 3-4 
transmission plots are laid out on the same set of hori- 
zontal axes. The dashed red lines denote T = 1, and 
solid curves are the transmission coefficient as defined in 
Eq. (23). The nonlinearity g increases in steps of 0.01 
in Fig. 4 as we move upward on the plot, and in steps 
of 0.1 in 3. Each transmission plot is at a convenient 
vertical offset for illustration, but the plots are not oth- 
erwise scaled or shifted. The nonlinearity is shown next 
to each curve for reference. As discussed in Section II C, 
the transmission coefficient may be greater than 1. This 
is not a novel physical feature; it arises due to the redef- 
inition of transmission as in Eq. (23) and the invalidity 
of superposition in this problem. 

Comparing the linear, i.e g — 0, transmission curve in 
Fig. 3 with the other transmission plots in the same fig- 
ure, we observe that there is a significant change in trans- 
mission behavior when we enter the nonlinear regime, 
even for very small nonlinearity, g ~ O(10~ 2 ). The most 
significant change occurs when /i is small, meaning that 
the potential barrier has a greater overall effect on the 
condensate. In the linear case, T oscillates about evenly 
on either side of unity for small /i, and the amplitude 
of oscillations is similar in either direction. In the non- 
linear case, we see that while transmission can still be 
less than unity, it does not drop as far below unity as in 
the linear case. We understand this difference to be due 
to the effective potential: for g = the operator in the 
NLSE is not biased, whereas for non-zero g it is. Thus, 



-|/Waaa/w^— 








10 



20 



30 



40 



FIG. 4: (Color online) Wide barrier, medium 
nonlinearity. Same as Fig. 3, but with steps of 
Ag = 0.1, from g = 0.1 to g = 1.0 An extended region 
of near-perfect invisibility of the barrier in the \x ~ 20 
to fi ~ 30 range translates to the right as the 
nonlinearity increases. 



by choosing parameters on the left hand side, we fix the 
effective potential and bias the system towards particular 
parameter sets. 

As g increases, new peaks appear in the small-^ regime, 
and the behavior of T changes significantly for smaller p. 
This regime does not appear for g > 0.2. The reason 
that the transmission curves depend strongly on g in the 
small- /i regime is that g and /j, are not completely in- 
dependent. For a given chemical potential p, when the 
nonlinearity g becomes larger than a certain cut-off value, 
solutions are generated that have p(x) £ C. This is in- 
valid because, when we solved the NLSE, we assumed 
p(x) £ M. We may take p(x) € M without loss of gener- 
ality, since Eq. (8) is a general polar representation of a 
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complex number. For large \x the cut-off in g is pushed 
to a region off the top of our plot; see figures in Sec. V 
for more details. The appearance of complex p(x) means 
that one or more of our assumptions are breaking down 
in this regime. A more detailed analysis of this situa- 
tion is considered in [37]. One such break-down is due 
to the appearance of quasi-bound states with complex 
H [45-47]. 

For larger values of n, the overall behavior of the trans- 
mission does not change significantly between each plot. 
However, we do see a shift in the transmission curve. 
For higher fj,, the transmission plot retains its shape and 
shifts to the right as g increases. This feature is espe- 
cially apparent when a computer is used to animate the 
transmission plots for increasing <?; in Fig. 4 we have 
endeavored to lay such an animation out on the page. 
In Fig. 4, the regime of the same curve-shape translat- 
ing gradually to the right with increasing g occurs when 

Another feature of note is that with the definition (23) 
of transmission, the amplitude of oscillations in T does 
not decrease monotonically as in the usual linear inter- 
pretation. In all of the transmission plots of Fig. 4, we 
see significant oscillations on either side of a transition 
region between [i ~ 20 and /i ~ 30. In this region, we 
have almost perfect resonance, i.e., T is very close to 
unity. Further analysis of transmission resonances will 
be presented in Sec. V. 



B. Narrow barrier 

T 

1-8 :vwvw>- ^^^ g_=_ 0.04 
1 .6 ; Vw ^^^-^ = _^^ B_= 0.03 

1.2 -vw^^^——^^- - 01 

i.o ^/\/\/\A^w_^-x^_fc= 

10 20 30 40 ^ 

FIG. 5: Transmission for narrow barrier with small 
nonlincarity. Note the difference in behavior between 
the linear and nonlinear regimes. 

To explore the narrow-barrier regime, we choose X2 — 
x\ = 0.1 and Vo = 10, so that the barrier is a factor 
of ten narrower than it is high; note that the height is 
the same as the wide-barrier case explored in Sec. IV A. 
Over the full range of values of fi and g that we consider, 
the narrow-barrier requirement, A 3> x x — x\, is satis- 
fied, as can be verified from Eq. (13); in the range we 
consider, A > 0.35. We first focus on transmission plots 



T 




5 10 15 20 25 30 35 



FIG. 6: Transmission for narrow barrier with increasing 
nonlinearity from bottom to top, from g = 0.1 to 
g = 2.0 in steps of Ag = 0.1. 



for small nonlinearity, in Fig. 5. The layout of the figure 
is the same as in Sec. IV A: five transmission plots are 
shown on the same set of horizontal axes, with vertical 
offsets in T for illustration, the plots are not otherwise 
shifted or scaled, and g is increased from zero in steps 
of Ag = 0.01. Comparing the linear g = transmission 
curve with those for g > 0, we again see a significant drop 
in the amplitude of oscillations as we go from the linear 
to the nonlinear regime. Again, the nonlinear transmis- 
sion curves tend to stay further above unity than below 
unity, although the difference is less pronounced than in 
the wide-barrier case. As a function of fi the transmission 
curve is much smoother as compared to Fig. 3. This is 
because there is always less than one wavelength fitting 
into the barrier, and thus a small change in /i cannot sud- 
denly cause an integer number of wavelengths to match 
the barrier width. 

In Fig. 6, we turn to the regime of medium nonlinearity, 
increasing g from 0.1 in steps of Ag = 0.1 up to g — 
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FIG. 7: Transmission resonances for narrow barrier. All 
behavior is described by regularly spaced straight lines 
except for very small /i, where some combinations of g 
and /i do not produce stationary solutions. 
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FIG. 8: Transmission resonances for wide barrier. 
Bifurcations occur in an extended region of near-perfect 
invisibility of the barrier. 



2.0. The amplitude of oscillations decreases smoothly 
for increasing fi, independent of g. There is no region 
of near-perfect resonance as in the wide-barrier case: T 
oscillates quasi-periodically about T = 1. In fact, the 
entire transmission curve simply translates smoothly to 
the right as g increases, with no real abrupt behavior even 
for the smallest values of /it allowed for each curve. Thus 
we can surmise that the abrupt behavior in Sec. IV A is 
due to the barrier width, and it is the width, not the 
height, that mainly controls the transmission curves. In 
the following section we will adduce further evidence to 
support this point. 



the resonances are sparse for low values of ii, with much 
more sporadic and abrupt behavior than the narrow bar- 
rier case. This is because an integer number of wave- 
lengths can fit into the barrier width, an effect which is 
more pronounced for small /z; for instance, for g = 0.02 
and n = 2, A = 2.27, while for fx = 40, A = 0.355. In 
the regions where lines of resonance appear, the spacing 
between these lines is smaller than in Fig. 7. The spac- 
ing increases as /i increases, though not as fast as in the 
narrow-barrier case. 



V. TRANSMISSION RESONANCES AND 
BIFURCATIONS 

To better understand the results of Sec. IV, in par- 
ticular the extended region of near-perfect invisibility of 
the barrier, we focus on the transmission resonances only, 
i.e., the points in the g—fi plane for which T = 1 to within 
a tolerance of 10~ 5 . 

We first treat the simpler case of the narrow barrier. 
The transmission resonances from Fig. 6 are shown in 
Fig. 7. This clearly shows that all transmission reso- 
nances translate smoothly to the right for increasing g, 
with a constant slope. This figure allows us to predict the 
location of all transmission resonances except for small 
/j,, where constraints in the allowed simultaneous values 
of g and /x cut off certain regions, as discussed in Sec. IV. 
We also observe in Fig. 7 that the spacing between curves 
increases with increasing /i. A functional fit to this spac- 
ing can allow us to predict transmission resonances in the 
entire g—fi plane. Since the spacing is the same starting 
from g = 0, it is straightforward to find such a function 
for any values of A, B, and 5q. 

With the narrow-barrier case as a reference, we move 
on to the more intricate wide-barrier case. Figure 8 shows 
the transmission resonances from Fig. 4. We observe that 



The resonances are sparse for lower values of fi and 
more uniform for higher values of /i. For mid-range fi, 
we see very different behavior. There are three regions 
between jj, — 20 and /i = 26 where the resonances are 
extremely dense. These correspond to the region of near- 
constant resonance seen in the transmission plots. This 
region shifts to the right as g increases, as observed in 
the transmission plots. Three bifurcations are present in 
this region. Bifurcations are typical of nonlinear systems; 
the particular example we show here for this parameter 
set is expected to be a generic feature for all parameter 
sets. Further numerical exploration found these three 
bifurcations continued in an apparently infinite sequence 
sloping away upwards to the right. We found that a linear 
function was not sufficient to fit such bifurcations, and 
they occurred near but not precisely at a value of A — 
1/2. We did not find bifurcations at other rational values 
of A, so this may be a coincidence; the barrier length in 
this parameter set is x-i — x\ — 10. We explored up to /x = 
70 and g — 10. Although our computationally intensive 
study was only performed thoroughly for this particular 
choice of parameters, we did spot checks through many 
regions of parameter space and observed similar extended 
regions of near-perfect invisibility of the barrier. 
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FIG. 9: Slopes of resonance transmission lines, as a func- 
tion of density offset, coarse sampling. The curves show 
B = 0.2, 0.5, 1.0, 1.5, 2.0 from bottom to top; a vertical 
offset is given to each curve for visualization, since they 
in fact all lie on top of each other. The points represent 
actual data, while the curves are a guide to the eye. 




1 ,,,, 2 ,,,, 3' ,,, 4 ,,,, 5 S 

FIG. 10: Slope of resonance curves as a function of input 
amplitude. Although the transmission resonances curves 
do not depend on B and 6o, they do depend strongly 
on A. Points are numerical calculations, while the solid 
curve is a least squares fit. 



VI. CONCLUSIONS 



We return to the narrow barrier case for further analy- 
sis. We consider the slope s of these parallel lines of reso- 
nance as a function of other physical parameters. For the 
parameter regimes considered in this study, the slope of 
the resonance lines does not depend on the value of the 
density offset, B, as seen in Fig. 9. In this plot, the slope 
is shown for several values of B. Each curve has been 
vertically shifted by a convenient offset for illustration; 
all curves in fact lie on top of each other. We note that 
the elliptic parameter m, which is strongly governed by 
the nonlincarity of the system, depends on both A and 
6, but does not depend on B, as described in Sec. II B. 
Similarly, we found that the slopes s do not depend on 
the spatial translational offset <5o- However, the slopes do 
depend strong on the amplitude A. The slope decreases 
as the amplitude of the input density increases. We find 
an exponential least-squares fit as shown in Fig. 10: 



s(Aj) = 0.167 + 0.485e-°- 549A ', (27) 



where we fixed Bj = 1.0. 



We developed a general method for obtaining station- 
ary states of the nonlinear Schrodinger equation for any 
piecewise constant potential. We applied this method to 
nonlinear scattering on a rectangular potential barrier, 
focusing on cnoidal waves, or dark soliton trains. This 
problem differs greatly from the textbook linear quantum 
mechanics problem of scattering on a rectangular barrier. 
Among novel nonlinear features are the following. First, 
the wavelength in the three regions (left of the barrier, 
on the barrier, and right of the barrier) need not be the 
same. As a consequence, such a barrier can be used to 
create one or more sharply localized dark solitons from a 
broad phonon-like input. Second, it is mainly the barrier 
width, not its area, that controls the kind of stationary 
states observed. For wide barriers in which the barrier 
width is larger than the wavelength extended regions of 
near-perfect transmission occur; the barrier is invisible 
over a range of interaction strength g and chemical po- 
tential fi. Third, in wide barriers an apparently infinite 
sequence of bifurcations appears in this invisibility re- 
gion. Despite the parameter space for solutions being 
large, we showed that it is just the amplitude that de- 
termines slopes of transmission resonance lines, greatly 
simplifying this complex problem. 

To relate our predictions to alternate physical input 
of average density, momentum, and energy of an atom 
laser formed from a Bose-Einstein condensate, the anal- 
ysis laid out in Ref. [12] may be used directly; thus we 
do not repeat it here. Future possibilities for this work 
include applying our method to various piecewise con- 
stant potentials of interest in applications, and a more 
detailed mathematical study of the bifurcations we found 
in transmission resonances. 

We acknowledge helpful discussions with Paul A. Mar- 
tin and Matthew Heller. This work was supported by the 
National Science Foundation under Grant PHY-0547845 
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Appendix A: Jacobi Elliptic Functions 

There are twelve Jacobi elliptic functions in all: sn, 
cn, dn, sc, nc, dc, cs, ds, ns [38]. Not all of these are in- 
dependent. The Jacobi functions are doubly-periodic in 
the complex plane and they depend on two parameters: 
an independent variable u and the elliptic parameter to. 
For real parameter m, we may assume < m < 1. If m 
is outside of this range, transformations can be made to 
write the Jacobi function in terms of functions whose pa- 
rameter is between and 1. For sn (u\m), which appears 
in the density p(x) of the BEC, these transformations 
are, for to < 0, 



sn(w|m) = 



1 — m 



1/2 



x sd 



1 — m 



, (Al) 



and for to > 1, 

sn (u\m) = TO~ 1/2 sn(wn 1/2 |TO _1 ). (A2) 

When in lies between zero and unity, the Jacobi func- 
tions can be interpreted geometrically as the analog of the 
hyperbolic and trigonometric functions. In this interpre- 
tation, the parameter m corresponds to the eccentricity 
of the ellipse. For m = 0, the Jacobi functions reduce to 
the trigonometric functions; for to = 1, to the hyperbolic 
functions. 

The Jacobi elliptic functions are defined as inverse in- 
tegrals [48] , and by the locations of zeros and poles in the 
complex plane [38]. They may be related to one another 
by various identities [49-51]. 



Proof. Compute the Wronskian of the two functions: 
W[sn p (u\m),sn q (u\m)] = 



sn p (w|m) sn 9 (w|TO) 
&K(u|m)] |;KHm)] 

gsn p (u|m,)sn 9-1 (u|m)cn (u|m)dn (u\m) 
n P-i 



(Bl) 



psn p (u\m) sn q (u\m) cn (u\m)dn (u\m) (B2) 

(B3) 



(q — p)cn (u\m)dn (u\m)sn p+q (u\m). 



The product of Jacobi elliptic functions is nonvanishing 
except on a set of measure zero; therefore, for p ^ q the 
functions are linearly independent. QED. 



2. Linear Independence of Products of Powers of 
sn (u\m) 



Let 



fi(x) = sn p (a|TO)sn 9 (w|m), 
/ 2 (x) = sn r (g | TO)sn s | to) , 



(B4) 
(B5) 



where p, q, r, s are integers, and a and u are linear func- 
tions of x. We may assume without loss of generality 
that a ^ u. The case a — u is analyzed in Section B 1. 
Computing the Wronskian of /i, / 2 , we find 

W(h,f 2 ) = 

(r — p)cn (a\m)dn (a|TO)sn r+p_1 (a|TO.)sn s+9 (ii|TO)+ 
(q — s)cn (u|m)dn (u|TO)sn s+9-1 (w|m)sn r+p (a|TO). (B6) 

The Wronskian vanishes Va, u only when r = p and q = s; 
that is, when /i and / 2 are not distinct functions. In all 
other cases, /i and / 2 are linearly independent. 



Appendix C: Numerical Convergence 

We consider several limits of the NLSE, Eq. (9). These 
are used to verify consistency of our code with known 
results. 



Appendix B: Proofs 

1. Linear Independence of Powers of sn (u\m) 

In the following we state a particularly vital proof 
which we have not found elsewhere in the literature, and 
is required to establish our exact solutions. Other proofs 
can found in Rcf. [37]. 

Theorem B.l. The functions sn p (u\m) and 
sn q (u\m) 7 p,q £ Z. are linearly independent for 



1. Linear Limit 

We consider the NLSE, Eq. (9), in the limit that g — > 
0. In this limit, the solution should reduce to the well- 
known linear scattering solution, which is analyzed in 
many elementary quantum mechanics texts. In the linear 
limit, we find that to — > 0, so that the linear density is 

p(x) = Asm 2 (bx + 6 ), (CI) 

since B = in the linear case [37]. For comparison with 
the nonlinear case, we define transmission as (pin) / (pi) . 
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The integral for (p) can be evaluated exactly in this case: 



(P) = 



2-rr/b 



dx A sin 2 (fa; + 5q) , 



\A. 



(C2) 
(C3) 
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Therefore, the transmission is 



T, = 



A 



in 



A, 



(C4) 



We can compare the value given by Eq. (C4) to the value 
T obtained numerically by the code. Transmission plots 
are shown in Figs. 11a and lib. A log plot of the error, 
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(c) 



FIG. 11: (a) Numerically computed transmission for lin- 
ear limit, (b) Exact linear value of transmission, (c) 
Error in transmission for linear limit. 



\T-T t \ 
T 



(C5) 



Constant Potential Limit 



is given in Fig. 11c. The maximum value of the error is 
0(1O -7 ), which is within the numerical tolerance of the 
code. 



Consider the NLSE, Eq. (9), with potential barrier 
(15), in the limit that Vo — > 0. In this case, boundary 
conditions are redundant and we expect all parameters 
to be constant Vx. By setting a "barrier" of Vo = in 
the code, we can verify that the code gives the correct 
solution; namely, the amplitude, period, and shifts in the 
density should not change at the "boundary" locations. 
Indeed this is the case, providing an additional verifica- 
tion of correctness for the code. A density plot for this 
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case, with a "barrier" of width 5 and height 1, is shown 
in Fig. 12. 

p(x) 
2.0h 



1.5 
1.0 
0.5 



-4 -2 2 4 6 8 10 
FIG. 12: Density plot for zero barrier. 

Since the density parameters do not change over space, 
we expect to find a transmission coefficient of 1. We 
compute and plot the error 



In 



\T-1\ 
T 

-I- ;iv<>- 



(C6) 



where T" aV g denotes the average value of transmission over 
the plot interval. A log plot of the error is given in Fig. 13. 
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FIG. 13: Error in transmission for zero barrier. 

Alternatively, we can consider a constant nonzero po- 
tential: V(x) = V c , Vx in Eq. (9). Again we expect 
all parameters to be constant Vx. We set a "barrier" of 
Vi = Vu = Viu = 2 and width 5 in the code and plot 
the density. The plot is shown in Fig. 14. 

Again, since the density parameters do not change over 
space, we expect to find a transmission coefficient of 1. 
We compute and plot the error 



In 



\T-l\ 
T 



(C7) 



where T aV g denotes the average value of transmission over 
the plot interval. A log plot of the error is given in Fig. 15. 

Therefore the code gives the expected results for den- 
sity and transmission in the constant-potential limits. 
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FIG. 14: Density plot for nonzero flat barrier. 
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FIG. 15: Error in transmission for nonzero flat barrier. 

3. Thomas-Fermi Limit 

We consider the limit that hd 2 ^ jdx 2 — > in the un- 
sealed NLSE, Eq. (1). In this limit, the unsealed NLSE 
becomes 

[g\*(x, t)\ 2 + V(x)] *(x, t) = ih^(x, t). (C8) 

Substituting Eq. (8) for the wave function ^ in (C8) and 
rearranging, we obtain 

gp 3 ' 2 + [V(x) - h^p 1 ' 2 = 0, (C9) 

so that either p(x) = 0, the trivial case, or else 

p(x) - [Hp - V(x)]/g, (C10) 

when g 0. Note that Eq. (C10) works for any spatially- 
dependent potential V(x). 

Equation (C10) is the well-known Thomas-Fermi 
limit [1] . It is relevant when the curvature of $ is nearly 
zero. This can be accomplished either by taking the limit 
A/B — >• 0, so that the amplitude of oscillations is small, 
or by taking the limit 1/6 —> oo, so that the wavelength 
of oscillations is large. 



14 



[1] F. Dalfovo, S. Giorgini, L. P. Pitaevskii, and S. Stringari, 
Rev. Mod. Phys. 71, 463 (1999). 

[2] W. Ketterle, D. S. Durfee, and D. M. Stamper-Kurn, in 
Proceedings of the International School of Physics "En- 
rico Fermi" (IOS Press, Amsterdam; Washington, D.C., 
1999), pp. 67-176. 

[3] B. P. Anderson and M. A. Kasevich, Science 282, 1686 
(1998). 

[4] E. W. Hagley, L. Deng, M. Kozuma, J. Wen, K. Hele- 
merson, S. L. Rolston, and W. D. Phillips, Science 283, 
1706 (1999). 

[5] I. Bloch, T. W. Hansen, and T. Esslinger, Phys. Rev. 

Lett. 82, 3008 (2000). 
[6] O. Zobay, E. V. Goldstein, and P. Meystre, Phys. Rev. 

A 60, 3999 (1999). 
[7] K. Helmerson, D. Hutchinson, K. Burnett, and W. D. 

Phillips, Phys. World 12, 31 (1999). 
[8] P. Engels and C. Atherton, Phys. Rev. Lett. 99, 160405 

(2007). 

[9] D. Dries, S. E. Pollack, J. M. Hitchcock, and R. G. Hulet, 

e-print arXiv: 1004. 1891 (2010). 
[10] V. Hakim, Phys. Rev. E 55, 2835 (1997). 
[11] N. Pavloff, Phys. Rev. A 66, 013610 (2002). 
[12] B. Seaman, L. D. Carr, and M. J. Holland, Phys. Rev. A 

71, 033609 (2005). 
[13] D. Witthaut, S. Mossmann, and H. J. Korsch, J. Phys. 

A: Math. Gen. 38, 1777 (2005). 
[14] B. Seaman, L. D. Carr, and M. J. Holland, Phys. Rev. A 

71, 033622 (2005). 

[15] T. Paul, P. Schlagheck, P. Leboeuf, and N. Pavloff, Phys. 
Rev. Lett. 98, 210602 (2007). 

[16] The results of our study can be qualitatively applied to a 
barrier of smooth shape with the same area as the rect- 
angular barrier. 

[17] M. Olshanii, Phys. Rev. Lett. 81, 938 (1998). 

[18] D. S. Petrov, G. V. Shlyapnikov, and J. T. M. Walraven, 
Phys. Rev. Lett. 85, 3745 (2000). 

[19] L. D. Carr, M. A. Leung, and W. P. Reinhardt, J. Phys. 
B: At. Mol. Opt. Phys. 33, 3983 (2000). 

[20] V. Dunjko, V. Lorent, and M. Olshanii, Phys. Rev. Lett. 
86, 5413 (2001). 

[21] Emergent Nonlinear Phenomena in Bose-Einstein Con- 
densates, edited by P. G. Kevrekidis, D. J. Frantzeskakis, 
and R. Carretero-Gonzalcz (Springer- Verlag, Berlin, 
2008). 

[22] L. D. Carr, C. W. Clark, and W. P. Reinhardt, Phys. 

Rev. A 62, 063610 (2000). 
[23] L. D. Carr, C. W. Clark, and W. P. Reinhardt, Phys. 

Rev. A 62, 063611 (2000). 
[24] B. Seaman, L. D. Carr, and M. J. Holland, Phys. Rev. A 

72, 033602 (2005). 

[25] L. D. Carr, K. Mahmud, and W. P. Reinhardt, Phys. 

Rev. A 64, 033603 (2001). 
[26] K. Rapedius, D. Witthaut, and H. J. Korsch, Phys. Rev. 

A 73, 033608 (2006). 
[27] H. A. Ishkhanyan and V. P. Krainov, Phys. Rev. A 80, 

045601 (2009). 

[28] T. Ernst and J. Brand, Phys. Rev. A 81, 033614 (2010). 
[29] C. Weiss and Y. Castin, Phys. Rev. Lett. 102, 010403 



(2009). 

[30] A. I. Streltsov, O. E. Alon, and L. S. Cederbaum, Phys. 
Rev. A 80, 043616 (2009). 

[31] J. Denschlag, J. E. Simsarian, D. L. Feder, C. W. Clark, 
L. A. Collins, J. Cubizolles, L. Deng, E. W. Hagley, K. 
Helmerson, W. P. Reinhardt, S. L. Rolston, B. I. Schnei- 
der, and W. D. Phillips, Science 287, 97 (2000). 

[32] A. Weller, J. P. Ronzheimer, C. Gross, J. Esteve, M. K. 
Oberthaler, D. J. Frantzeskakis, G. Theocharis, and P. G. 
Kevrekidis, Phys. Rev. Lett. 101, 130401 (2008). 

[33] S. Burger, K. Bongs, S. Dettmer, W. Ertmer, K. Sen- 
gstock, A. Sanpera, G. V. Shlyapnikov, and M. Lewen- 
stein, Phys. Rev. Lett. 83, 5198 (1999). 

[34] L. Khaykovich, F. Schrcck, F. Ferrari, T. Bourdel, J. Cu- 
bizolles, L. D. Carr, Y. Castin, and C. Salomon, Science 
296, 1290 (2002). 

[35] W. P. Reinhardt and C. W. Clark, J. Phys. B: At. Mol. 
Opt. Phys. 30, L785 (1997). 

[36] S. Stellmer, C. Becker, P. Soltan-Panahi, E.-M. Richter, 
S. Dorscher, M. Baumert, J. Kronjager, K. Bongs, and 
K. Sengstock, Phys. Rev. Lett. 101, 120406 (2008). 

[37] R. R. Miller, Master's thesis, Colorado School of Mines, 
2010. 

[38] Handbook of Mathematical Functions, edited by M. 
Abramowitz and I. A. Stegun (National Bureau of Stan- 
dards, Washington, D. C, 1964). 

[39] The final condition, /i~ = is not strictly required 
when the density vanishes on the boundary; however, in 
the latter case, one can argue this condition holds on 
physical grounds [37]. 

[40] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. 
Flannery, Numerical Recipes in C: The Art of Scientific 
Computing (Cambridge Univ. Press, Cambridge, U.K., 
1993). 

[41] K. E. Strecker, G. B. Partridge, A. G. Truscott, and R. G. 

Hulet, Nature 417, 150 (2002). 
[42] B. A. Kalinikos, M. M. Scott, and C. E. Patton, Phys. 

Rev. Lett. 84, 4697 (2000). 
[43] M. Wu, B. A. Kalinikos, and C. E. Patton, Phys. Rev. 

Lett. 93, 157207 (2004). 
[44] Z. Wang, A. Hagerstrom, J. Q. Anderson, W. Tong, M. 

Wu, L. D. Carr, R. Eykholt, , and B. Kalinikos, Phys. 

Rev. Lett, to be submitted (2010). 
[45] N. Moiseyev, L. D. Carr, B. A. Malomed, and Y. B. Band, 

J. Phys. B: At. Mol. Opt. Phys. 37, LI (2004). 
[46] L. D. Carr, M. J. Holland, and B. A. Malomed, J. Phys. 

B: At. Mol. Opt. 38, 3217 (2005). 
[47] G. Dekel, V. Farberovich, V. Fleurov, and A. Soffer, 

arXiv e-print:0911.1537 (2009). 
[48] F. Bowman, Introduction to Elliptic Functions, with Ap- 
plications (Dover, New York, 1961). 
[49] H. E. Fettis, Mathematics of Computation 26, 965 

(1972). 

[50] A. Khare, A. Lakshminarayan, and U. Sukhatme, Pra- 
mana Journal of Physics 62, 1201 (2004). 

[51] B. Dayton, in Theory of Equations (Oakton, Des Plaines, 
IL, 2002), Chap. Analysis - Elliptic Functions. 



